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We present a theoretical study on the de Haas-van Alphen (dHvA) oscillation in the vortex state 
of type-II superconductors, with a special focus on the connection between the gap anisotropy and 
the oscillation damping. Numerical calculations for three different gap structures clearly indicate 
that the average gap along the extremal orbit is relevant for the magnitude of the extra damping, 
thereby providing support for experimental efforts to probe gap anisotropy through the dHvA signal. 
We also derive an analytic formula for the extra damping which will be useful to estimate angle- 
and/or band-dependent gap amplitudes. 
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A considerable number of materials have been found 
in the 1990s to exhibit de Haas-van Alphen (dHvA) os- 
cillation in the superconducting vortex state || , the phe- 
nomenon discovered by Graebner and Robbins in 2H- 
NbSe2 ||. Many theories have been put forward during 
the same period to explain this fundamental phenomenon 
observed in the system without a well-defined Fermi sur- 
face Q. However, there is yet no established theory com- 
parable with the normal-state Lifshitz-Kosevich (LK) 
theory ||. Moreover, what is lacking seems to be a clear 
physical picture of the mechanism of the extra oscillation 
damping in the vortex state. To improve the situation, 
we present here a numerical study combined with an an- 
alytical one. 

Novel aspects in our study are summarized as follows: 

(i) We perform three-dimensional numerical calcula- 
tions for the dHvA oscillations in the vortex state us- 
ing the Bogoliubov-de Gennes equations. To date, nu- 
merical studies have been carried out only for the two- 
dimensional s-wave model, such as that of Norman et al. 
(NMA) 0. Our results show that dHvA oscillations in 
three dimensions are quantitatively different in the vor- 
tex state from those in two dimensions. 

(ii) To clarify the connection between the oscillation at- 
tenuation and the gap anisotropy, we use three different 
gap structures: s-wave, <i-wave with four point nodes in 
the extremal orbit, and p-wave with a line node in the ex- 
tremal orbit; see Eq. below, ft is thereby shown that 
the attenuation is determined by the average gap along 
the extremal orbit in zero magnetic field (see Fig. 1), in 
disagreement with the theory of Miyake || . This result 
still indicates that one can probe the angle- and/or band- 
dependent gap amplitude in zero field with the dHvA ef- 
fect in the vortex state. The origin of the discrepancy 
from Miyake's result is discussed in some detail. 

(iii) We derive a new analytic formula for the extra 
Dingle temperature in the vortex state: 



k B T A =0.5f(|A p | 



m b c 1-B/H c2 
' neh B 



(1) 



(|A p | 2 ) co denotes the average gap along the extremal or- 
bit in zero field, is the band mass, and B and H C 2 
are the average flux density and the upper critical field, 
respectively. 

Equation ([!]) is derived below through the second-order 
perturbation with respect to the pair potential, i.e., an 
approach from -ff c2 , which is useful to estimate the gap 
amplitude along the extremal orbit. A difference from 
Maki's formula || lies in the prefactor where the Fermi 
velocity v-p is absent. Indeed, a dimensional analysis on 
the second-order perturbation tells us that the Landau- 
level broadening in the vortex state should be of order 
\M-°\B)\ 2 /Tilo c , where huj c is the cyclotron energy and 
&°\B) cx y/(\A p \ 2 ) co (l-B/H c2 ) is essentially the aver- 
age gap along the extremal orbit. This leads to Eq. (|l|) 
except for the numerical constant. 

Terashima et al. reported a dHvA experiment on 
YM2B2C where an oscillation (labeled a) is seen to per- 
sist down to a field of order 0.2if c2 . On the other hand, 
a specific-heat experiment at H = shows a power-law 
behavior oc T 3 at low temperatures S , indicating the ex- 
istence of the gap anisotropy in this material. This gap 
anisotropy may also play an important role for the dHvA 
signal far below iJ c2 . By using Eq. (|l|), the average gap 
of this band a will be estimated at the end of the paper. 

Model: Our starting point is the Bogoliubov-de 
Gennes equation for the quasiparticle wavefunctions u s 
and v* labeled by a quantum number s with a positive 
eigenvalue E s : 



dr 2 



n{r u r 2 ) A(n,r 2 ) 
-A*(ri,r 2 ) -W*(ri,r 2 ) 

u s (ri) 

-<( r i) 



u 5 (r 2 ) 
-v*(r 2 ) 



(2) 



Here, T = 0.125 is a dimensionless quantity characterizing 
the Landau-level broadening due to the pair potential, 



Here A is the pair potential and 7£ denotes the normal- 
state Hamiltonian in the external field H |j z; both are 
2x2 matrices to describe the spin degrees of freedom. 
We adopt as H the free-particle Hamiltonian. As for the 
pair potential, we wish to consider three cases to yield 
the following energy gaps at H = Q: 
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where p is the momentum, W p denotes some cut-off func- 
tion with W PF = 1 (pp: Fermi momentum), 9 p (ip p ) is 
the polar (azimuthal) angle, and o^-'s are the Pauli ma- 
trices. The first one is the isotropic s-wave state, whereas 
the latter two have four point nodes (d-wave) and a line 
node (p-wave), respectively, in the extremal orbit of the 
xy plane perpendicular to H || z. It then turns out that 
the orbital part of the corresponding pair potentials in 
finite fields can be expanded with respect to r = i 1 ! — r 2 
and R=i(r!+r 2 ) as [§ 



A(r 1; r 2 ) =4E A(JVc) Wg q (Rx) £ (-1) 



V2 



N T mp z 



W p cos p 
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Here N c and N r denote the Landau levels in the average 
flux density J5, q is an arbitrarily chosen center-of-mass 
magnetic Bloch vector, and m signifies the relative angu- 
lar momentum along the z axis so that m = 0, 0, and ±2 
for the three cases, respectively, in a system with M 2 /2 
flux quanta and the length L z along the z axis. The ar- 
guments rj_ and Rj^ denote the xy components, and p 
and 9 P are to be evaluated at p= (p 2 z + h 2 N r /l"^) 1 / 2 and 

6» p = tan- 1 h ^J l z with iB-ihc/eB) 1 / 2 . See ref. § for 

the expressions of the basis functions and ipff i7n . 

A significant advantage of Eq. (Q) is that the coeffi- 
cients {A^ Ar<: ' ) }^ =0 completely specify the pair poten- 
tial, and the first few terms suffice to describe those of 
BZ0.1H c2 The self-consistent mean-field theory 

P, fl^|]ll| ] predicts an oscillatory reentrant behavior of H c2 
and A^ c ). On the contrary, such a singular behavior of 
H c2 has never been identified definitely in any materi- 
als displaying the dHvA oscillation in the vortex state. 
Leaving the puzzling discrepancy for a future study, we 
here adopt the quasiclassical £>S N ^ rather than the fully 
self-consistent one, with the square-root behavior 



m=a{l-B/H c2 



,1/2 
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of the mean-field second-order transition for the domi- 
nant N c = level. Then the best choice for {A'^)}^ _ 
would be to use the results from the Eilenberger equa- 
tions. Since our main interest lies in studying the dif- 
ferences in the oscillation damping among various gap 
structures, however, we here adopt a model form of 
^( N o) determined by requiring that the maximum of 
^JdR \Jdr A(n, v 2 )e~ ip - r,n \ 2 be equal to h? Q {l-B/H c2 ), 
where V is the volume of the system and Ao denotes 
the maximum energy gap of the weak-coupling theory at 



T = H = 0. This is possible within the lowest-Landau- 
level approximation of retaining only which is ex- 
cellent for B^0.lH c2 The resulting A<°) also 
displays the square-root behavior, and our numerical cal- 
culation shows that a 2 « 0.5 A 2 in all the three cases. We 
are planning to report on the best choice for a from the 
Eilenberger equations. 

The use of the quasiclassical pair potential (||) has an 
advantage: Self-consistent calculations for the continuum 
model y necessarily result in a rather small number of 
Landau levels below the Fermi energy Ep, i.e., Np ~ 10 
at H c2 , hence failing to meet the condition Np ^> 1 ap- 
propriate for real materials. The present calculation with 
Eq. (|^) is free from such a limitation. We have performed 
calculations including about 50 Landau levels below ep 
for the extremal orbit at H c2 . 

Since the relevant materials have large Ginzburg- 
Landau parameter k>1, we also neglect the screening 
in the magnetic field. Indeed, the effect has been shown 
to be irrelevant for the oscillation damping We put 
T = 0, adopt Wp = e-(^/°- l£F ) 4 with £ p the normal-state 
one-particle energy measured from sp, and choose the 
cyclotron energy at H c2 as hu> c2 = k^T c , in accordance 
with hu} c2 /k^T r = \^Z for real materials. 

Numerical Method: To solve it numerically for 
the above model, we transform Eq. (^) into the eigen- 
value problem for the expansion coefficients of u s and 
v s in the quasiparticle basis functions {ij^Nka}, where 
k is a quasiparticle magnetic Bloch vector and a (— 
1,2) signifies twofold degeneracy of the orbital states 
Then it can be solved separately for each ka 
due to the translational symmetry of the vortex lattice. 
The overlap integral between i/ , ^oq(-^' J -)^w!m( rj -) anc ^ 
ipN 1 k iai {ri±)^N2k 2 a 2 ( r 2±) vanishes unless q = k x + k 2 
and a.\ — a 2 so that the calculations are simplified 
greatly. The corresponding eigenstate is labeled explic- 
itly by s= (v\iap z <7) with v (a) the band (spin) index. 

Numerical results: Figure 1(a) presents the oscilla- 
tion of the s-wave magnetization as compared with the 
normal-state one. We see clearly that the oscillation 
frequency is unchanged from that of the normal state. 
With haj c = k^Tc at H c2 , the oscillation is observed to 
persist down to a rather low field of H c2 / B SSl.8, i.e., 
B^,0.55H c2 , which is lower than 0.8H c2 where fiuj c be- 
comes equal to the spatial average of the energy gap: 
Ao(l — B/H^) 1 / 2 . This is partly because the gap is 
smaller within the extremal orbit, as shown by Brandt et 
at [ fl3| The points with error bars in Fig. 1(b) comprise 
the corresponding Dingle plot for the extra damping fac- 
tor R s obtained by numerical differentiation. This extra 
damping at high fields shows the behavior cx l—B/H c2 in 
the logarithmic scale, but an irregularity sets in around 
0.55-ff c2 where the oscillation disappears. We attribute 
this irregularity to the effect of the bound-state forma- 
tion in the core region. The lines are the predictions from 
various theoretical formulas. The Maki formula M repro- 
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FIG. 1. Left figures: the oscillatory part M osc of the mag- 
netization in the vortex state (full lines) as compared with 
the normal-state one (dotted lines) for the three cases of Eq. 
(^). Right figures: the corresponding Dingle plots (points 
with error bars) compared with theoretical predictions. 

duces the correct functional behavior oc 1—B/H C 2 at high 
fields, but the prefactor is too large. The NMA formula 
j| , deduced from the two-dimensional self-consistent nu- 
merical results with Np ~ 10 at H C 2, clearly shows a more 
rapid attenuation incompatible with our numerical data. 
The theory of Dukan and Tesanovic |14|, which would 
predict R s = in the clean limit of T = 0, is also incon- 
sistent with the data. The solid line is due to Eq. (|l|) 
derived below; it reproduces both the functional depen- 
dence and the attenuation magnitude of the numerical 
data. 

The above discrepancy between our numerical results 
and those of NMA may be attributed mainly to a dif- 
ference in dimension. Indeed, the dHvA signal in three 
dimensions differs from that of two dimensions in that 
some finite region Sp z around the extremal orbit is rel- 
evant. Most of the Landau levels in the region do not 
satisfy the particle- hole symmetry with respect to ep, so 
that the effect of the pair potential is rather small to be 
handled with the second-order perturbation. This may 
be the reason why the oscillation is more persistent in 
our three-dimensional numerical results. This point is 
crucial in our derivation of the analytic formula below. 



Figures 1(c) and 1(d) show the results for the second 
energy gap in Eq. (||) which vanishes at four points on the 
Fermi surface along the extremal orbit. The damping is 
seen to be strong and not much different from the s-wave 
case. Equation (fil) gives a good fit to the numerical data 
for H C 2/B£ 1.8, thereby indicating that the average gap 
along the extremal orbit is relevant for the extra damp- 
ing. This result is in disagreement with Miyake's theory 
that point nodes in the extremal orbit should weaken the 
attenuation |j| . Indeed, Miyake obtained his analytic for- 
mula by applying a semiclassical quantization condition 
for the expression of the electron number N at B = 0. 
Neither of his starting point N(B = 0) nor the use of the 
quantization condition may be justified for describing the 
dHvA signals observed mainly near H C 2- 

In the line- node case of Figs. 1(e) and 1(f), in con- 
trast, we observe a much weaker damping in favor of 
Miyake's idea j5|. However, the non-zero extra damp- 
ing can only be explained by considering the contribu- 
tion of some finite width near the extremal orbit. Using 
Eq. (|l|) to obtain the best fit to the numerical data, i.e., 
the solid line of Fig. 1(f), we estimate that the region 
|Pz| SSpf sin^ contributes to the extra attenuation. This 
is in rough agreement with the estimation from the Fres- 
nel integral J"_^exp[— i(^2ttNf p z /pf) 2 ] dp z which ap- 
pears in obtaining the LK formula: Cutting the infinite 
integral at \/2ttNf\p z \/pf ~ 10 with Nf ~ 50 yields a 
similar value for the relevant range of p z . 

Analytic Formula: The Luttinger-Ward free-energy 
functional corresponding to Eq. (0) is given by pq| 



= 



knT 



Tr In 



H-ienI A 
-A* -H*-ie n l 

e ie "°+l 
e - 4£ "°+l 



+ •••, (6) 



where e n /h is the Matsubara frequency and 0+ is an 
infinitesimal positive constant. The terms • • • may be 
expressed solely with respect to the pair potential so that 
they can be neglected in the present model to consider 
the oscillatory part. This may be transformed into 



s 



k B Tln(l+e- E ^ k * T ) + E S J |v s (r)| 2 dr 



(7) 



Since we are interested in the extra damping in the vortex 
state, we adopt as E s and v s the expressions from the 
second-order perturbation with respect to A. They are 

ENkap z a = l&Vp^l+^kp/ign^jVp^) , (8) 

JW N k ap ,a(r)\ 2 dr= 0(-(, Npz(T ) + i 1 ( fil pz sign(£, Npzr7 ) , (9) 
where 6>(£) is the step function and = J2n'\II 

V^ka( r l^)^q-ka( r 2^) C "" < r^ 2V ' i A ( r l: r 2)rfr 1 dr 2 | 2 / 

(6vp z er + 6v-p 2 -o')" • The first terms on the right-hand 
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side of these equations are just the normal-state expres- 
sions. The second terms, on the other hand, denote the 
finite quasiparticle dispersion in the magnetic Brillouin 
zone and the smearing of the Fermi surface, respectively, 
due to the scattering by the growing pair potential. It 

is useful to express f]NUp z m terms of A( \B) and the 
cyclotron energy fiu) c of the extremal orbit as 



The quantity ^/J^L thus defined is dimensionless, and 
we realize that the main B dependence in Eq. (|l^) lies in 
the prefactor | U°\B)\ 2 / {hu c ) n . The explicit expression 

of vnI Pz is s iven b y 



\&°XB)\ Q 
(hu> c ) n 



An) 



(10) 



4 ^ 

N'mm' 



_(„) _ Aff \(NN'\0N + N')\ 2 (N+N'+m\2k-q) 

VNkp » ~ ~ [N+N'-2{N F + S)] n 

x{2k-q\N+N'+m') x { sin 4 6» p , (11) 

cos 2 p 

where the overlap integrals are given by eqs. (3.23) and 
(3.29) of ref. §, S = S{B, Pz ) (\S\ < 1/2) specifies the lo- 
cation of sf between the two closest Landau levels, and 
m, m! = 0, 0, and ±2 for the three cases of Eq. ([I]), respec- 
tively. The corresponding normalized density of states: 



D 



(") r~\ 



* ka 



~(n) 



), 



(12) 



will play a central role in the following. 

Substituting eqs. (||) and (pi) into Eq. (fit), we find that 

(2) 

the terms containing r] N ^ p may be neglected due to the 
cancellation between the particle and hole contributions. 
The remaining term can be transformed with the stan- 
dard procedure. We thereby obtain, for the first har- 
monic of Q/V, the expression: 



ill k^T ^ 



dNcos(2irN) / dp 



-1/2 



dfj 



(13) 



The function D K N ' {fj) depends on (N,p z ), but may be 

replaced by a representative one of (fj) to be placed 
outside the N and p z integrals fHI , where the recov- 
ered index I specifies the s-, d-, or p-wave case of Eq. 
(|bT|). It may also be acceptable to use a Lorenzian for 

it: D { l\fi) = fi/ir{fj 2 + T 2 e ) @. We thereby obtain an 
expression for the magnetization which carries an extra 
damping factor: 



Rs(B): 



D^\fj) exp 



-27rt<7|AW(B)| a /(fiWc) 5 



= exp[-27rf,|A(°)(B)| 2 /(^c) 2 



dfj 
(14) 



Thus, the superconductivity gives rise to an extra Dingle 
temperature of k B T A = f £ \M°\B)\ 2 /tt?iuj c , or equiva- 
lently, the extra scattering rate of t^ 1 = 2irk-QT A /h. 

Equation (^) may have an advantage that one can 
trace the origin of the extra dHvA damping definitely 
to the growing pair potential, which brings about a fi- 
nite quasiparticle dispersion ( |Tl| ) in the magnetic Bril- 
louin zone and the corresponding Landau-level broaden- 
ing, Eq. (|l^). Moreover, Eq. ( pi] ) reveals that this broad- 
ening near H C 2 is closely connected with the zero- field gap 
structure, Eq. (J3|) . 

There seems to be no analytic way to estimate T s , 
so we determine it using the best fit to the numerical 
data of Fig. 1(b). Using Eq. (§) with a 2 = 0.5A 2 ,, the 
procedure yields T = T s = 0.125, as noted before. It is 
also clear for the anisotropic cases that the average gap 
around the extremal orbit is relevant for the extra atten- 
uation, as may be realized from Eq. (|TT|) . We hence put 



a 2 F^0.5(|A p | 2 ) oo r s . We thereby obtain Eq. (g), which 
gives a good fit to the d-wave numerical data without any 
adjustable parameters; see Fig. 1(d). 

Concluding Remarks: We have shown explicitly 
that the dHvA effect in the vortex state can be a pow- 
erful tool to probe the average gap along the extremal 
orbit. Such an experiment has recently been performed 
on UPd 2 Al 3 by Inada et al. |J], and Eq. (Q) will be 
useful in similar experiments to estimate angle- and/or 
band-dependent gap amplitudes. Using the a oscillation 
of YM2B2C observed by Terashima et al. @, for exam- 
ple, we obtain (|A p | 2 ) 0O = 1.5meV for this a band, a 
value much smaller than 2.5meV from the specific-heat 
measurement ||. It should be noted, however, that the 
numerical factor 0.5 in Eq. (Q) comes from our model 
described below Eq. (||) . It has to be replaced by the re- 
sult from the Eilenberger equations. We are planning to 
report on this together with detailed comparisons with 
experiments in the near future. 
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